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We present three methods for calculating the T-matrix for modeling arbitrarily shaped micro- 
sized objects. When applied to microrotors, their rotation and mirror symmetry can be exploited to 
reduce memory usage and calculation time by 2 orders of magnitude. In the methods where the T- 
matrix elements are calculated using point matching with vector spherical fields, mode redundancy 
can be exploited to reduce calculation time. 

PACS numbers: Valid PACS appear here 



I. INTRODUCTION 

• Motivation: optical rotors 

• DDA [? ? ] 

• Memory requiements of the A matrix, swapping 

• Compressed A matrix 

• T-matrix [? ] 

• DDA T-matrix [? ] 
Our T-matrix methods are 

• near VSWF field point matching with DDA field 

• far VSWF field point matching with DDA field 

• rotate [? ] and translate [? ] 



II. OPTIMIZING THE DISCRETE DIPOLE 
APPROXIMATION INTERACTION MATRIX 

The size of the microdevices we model may exceed 10- 
20 wavelengths in size, which may well require compu- 
tational time in excess of several days and RAM be- 
yond that available. To circumvent these limitations, 
we exploit the discrete rotational and/or mirror symme- 
try of a microcomponent. This is closely tied with the 
link between DDA and the T-matrix method. In the 
T-matrix method, the fields are represented as sums of 
vector spherical wavefunctions (VSWFs) [? ? ], and to 
use DDA to calculate a T-matrix, we can simply calcu- 
late the scattered field (and its VSWF representation) for 
each possible incident single-mode VSWF field in turn. 
The important point is that each VSWF is characterized 
by a simple azimuthal dependence of exp(im(f>), where m 
is the azimuthal mode index. If we consider a group of 
dipoles that are rotationally symmetric about the vertical 



FIG. 1: 



axis — in the case of figure 1, there is 4th-order rotational 
symmetry — the magnitude of the incident field will be 
the same, the field differing by the vector rotation 6 and 
phase factor of exp(irn0); only the dipole moment of one 
repeating unit of the total number of dipoles needs to be 
known. In spherical coordinates, 



sph psph 



P ± p exp(im(j) q 



(i) 



However, our implementation of the interaction matrix[? 
], electric field and dipole moments were in cartesian co- 
ordinates system(what was the advantage?). This re- 
quires a transformation of the dipole moment at the first 
segment from the cartesian to spherical coordinate sys- 
tem followed by a transformation back from spherical to 
cartesian at the rotational counterpart dipole, 



P c q art = C q SP? art exp(im<f> q ), 
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is an orthogonal matrix that transforms the P£ art vector 
into the spherical coordinate system and 
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parts: 



FIG. 2: (a) 8-dipole example, (b) The 2 dipoles required to 
completely specify all dipole moments. 



is a transpose of S but at the azimuthal coordinate of the 
rotational counterpart dipole. 

This brings up question of how to reduce the num- 
ber of equations such that only one rotational unit needs 
to be solved. Conventionally, the interaction matrix as 
defined in [? ] represents the coupling between each 
dipole with all other dipoles. Figure 3(a) shows the in- 
teraction matrix for the example set of dipoles shown in 
figure 2(a). In general the matrix will be made up of 
N x N cells for N dipoles; each cell is a 3 x 3 tensor. In 
the example, the matrix is made up of 8 x 8 cells. A diag- 
onal cell represents the self interaction (or the inverse of 
the polarizability) and an off-diagonal cell represents the 
coupling between different dipoles. Taking advantage of 
the equal amplitudes and known phase factors between a 
dipole and its rotational counterparts, we can reduce the 
interaction matrix. Taking the example in figure 2(b), 
we contruct the interaction matrix as if there were only 
2 dipoles but we aggregate the contribution from the ap- 
propriate dipoles. For the off-diagonal cells, the coupling 
between a dipole with the other dipoles including their 
rotational counterparts are summed as follows 
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where Q is the order of discrete rotational symmetry, 
m is the azimuthal mode of the incident VSWF field, q 
is the rotational segment number, the rotational angle 



2nq/Q, 



M 
r jk 



is the distance from points to the 



rotationally symmetric points 
vector from points rj to rl . 
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and is the unit 



The coordinate for a given 
rotational symmetric point is calculated using a rotation 
about the z-axis, 
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For the diagonal cells, the "self interaction" includes 
the coupling between a dipole and its rotational counter- 
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Figure 3(b) shows the interaction matrix representation 
for the example dipole system in figure 2(c). The com- 
pressed interaction matrix is a factor of Q 2 smaller than 
the conventional matrix. Having precalculated the inci- 
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FIG. 3: (a) Full interaction matrix, (b) Symmetry reduced 
interaction matrix 

dent fields Ej. inc at each dipole of the rotational unit, 
we solve for the dipole moments Pj for the dipoles with 
a reduced set of linear equations. The dipole moments 
and fields of the rotational counterpart dipoles can be 
calculated by applying the rotational matrix 6 and phase 
factor exp(im(j)). 

We can exploit mirror symmetry in a similar fashion, 
since the VSWFs possess either even or odd parity w.r.t. 
the xy-plane. This allows a reduction in size by a further 
factor of 4. 



III. METHODS FOR CALCULATING THE 
T-MATRIX 

A. Near field point matching 

• vswf field matched with DDA field 

• solve for p's and q's 

• contract T-matrix column by column 

Pnm ,TEnm, (8) 

Qnm — ,TMnm (9) 



FIG. 4: (a) Near field matching, (b) Octant near field match- 
ing. 
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B. Far field matching 

• DDA farficld 

• VSH farficld 

• solving for p nm and q nm 

• Constructing the T-matrix 



Far field matching 



M^Hkr) = ^(Ti) n+1 exp(±ikr)C nm (9^) (10) 
N^\kr) = ^( T ir +1 exp(±ikr)B nm (0 }( t>) (11) 
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e) 







FIG. 5: a) 2 = translation b) z onto y axis rotation c) z' = 
translation d) z' onto a; axis rotation e) 2" = translation 
f) rotations back onto the original z-axis g) original orienta- 
tion 



IV. MODE REDUNDANCY AND THE ORDER 
OF DISCRETE ROTATIONAL SYMMETRY 

• Motivation, computational savings 

• Explain mth order discrete rotational symmetry 
and VSWF m modes 

• redundant modes 

• Floquet's theorem 

• mscat = mine +ip (5) 

• For a cube, p — 4, and mscat = mine, mine 4, mine 



Pnm 
Qnm 



Q>nm 



(15) 
(16) 



C. Rotation and translation of vector fields 

• Videen translations [? ] , precalculated 

• Rotations of axes 

• Constructing the T-matrix, cycle through m's and 
n's 



• time saving 



RESULTS 



A. Phase functions for the spheres and cubes 

comparision 

• Sphere Mie soln 

• cube point matching [? ] (toolbox) 

• also Wriedt f? 1 
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FIG. 6: The phase function for cubes with the refractive index 
of 1.5 and widths of a) 0.75A b) 1.25A. 



B. Cross rotor torque calculations 

• convergence of torque 

• low nrel so A/5 suffices 

• reference Theo's results 

VI. CONCLUSION 

Combination of optimization techniques 

• DDA rotational and mirror symm, interaction ma- 
trix compression 

• Far field and near field matrix, octant matching 
grid points 

• exploiting mode redundancy 
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Abstract 

We present a method of incorporating the discrete dipole approximation 
(DDA) method with the point matching method to formulate the T-matrix 
for modeling arbitrarily shaped micro-sized objects. The T-matrix elements 
are calculated using point matching between fields calculated using vector 
spherical wave functions and DDA. When applied to microrotors, their dis- 
crete rotational and mirror symmetries can be exploited to reduce memory 
usage and calculation time by orders of magnitude; a number of optimiza- 
tion methods can be employed based on the knowledge of the relationship 
between the azimuthal mode and phase at each discrete rotational point, and 
mode redundancy from Floquet's theorem. A 'reduced-mode' T-matrix can 
also be calculated if the illumination conditions are known. 

Key words: DDA, T-matrix, discrete rotational symmetry, point matching 



1. Introduction 

Optical tweezers [jj can be used to exert forces and torques and thus drive 
micromachines 0, 0] • This opens up a new field of micro engineering, whose 
potential has yet to be fully realized. To aid in designing micromachines, 
we employ a number of modelling methods ranging from the generalized 
Lorenz-Mie theory (GLMT) 3, the FDFD/ T-matrix hydrid method Q to 
the discrete dipole approximation (DDA) [y, [7|. We use the abovementioned 
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methods to formulate the T-matrix because repeated calculations with differ- 
ent illumination conditions are required. The extended boundary condition 
method (EBCM) @, [9[ is most commonly used to calculate the T-matrix; 
however, numerical problems are encountered with high aspect ratio struc- 
tures [lOt ] . 



Given the structure of typical optical micromachines [111 ], we have found 



the DDA method most suitable. In DDA, the scattering object is repre- 
sented as a collection of dipole scatterers, and the total scattering problem, 
including the coupling between the dipoles, is solved. DDA is well-suited to 
modelling optical micromachines. Firstly, only the volume of the actual par- 
ticle needs to be discretized, while both the particle and surrounding medium 
in a volume enclosing the particle are discretized in other general methods 
such as the finite-difference time-domain method (FDTD) and finite element 
methods (FEM). Considering that structures as shown in figured^ are not 
unusual, where the particle occupies only a relatively small fraction of the 
nearby volume, this can mean a considerable saving in required memory and 
time. Secondly, DDA performs well for relatively low contrast scatterers, 
which is typical of most optical micromachines so far, usually constructed 
from a polymer material 12j and deployed in a dielectric liquid. Thirdly, it 
is relatively simple to obtain the T-matrix via DDA if repeated calculations 
are desired 131 ] . Finally, it is possible to exploit discrete rotational symme- 



try of a particle to reduce the computational resources, including both time 
and memory, by orders of magnitude. This last factor is important, since 
optical micromachines are often large in overall dimension compared to the 
wavelength (while having wavelength scale features forcing the use of electro- 
magnetic theory rather than geometric optics), and the available resources 



can place such devices beyond practical calculation [14]. Figure [T] shows a 
typical case. 

By itself, the DDA method does not give us the T-matrix. One method of 
formulating the T-matrix via DDA 13| is to transform the field contributions 
of each dipole to aggregate their contributions at a common origin. We 
embarked on a different method — incorporated the DDA with the point 
matching method; both the near and far-field point matching methods are 
discussed in this paper but only the results of the former are presented. We 
will also discuss the optimization methods employed. 
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(a) (b) (c) 

Figure 1: (a) An optically-driven microrotor. (b) Discrete dipole representation of the 
rotor, (c) Discrete rotational and mirror symmetries allow the modelling of only the 
repeated segment, with a major reduction in required memory and time. 



2. DDA revisited 



With DDA, the microrotor in figure UK can be represented by a set of 
dipoles in a cubic lattice, as in figure [lb, fashioned to approximate the tar- 
get structure. The lattice spacing has to be sufficiently small, not only to 
accurately represent the structure — stair-casing is becomes a problem for 
curved surfaces — but small enough compared to the wavelength of the in- 
cident light; the criterion is specified in |7|. 

Ei nCi j, at each dipole and the 



We first have to calculate the incident field, 
interaction matrix [?J 
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jk 



exp(ikr jk ) 



k 2 (f jk f jk - 1 3 ) + l ^Jt — -(3r jk r jk - 1 3 ) 
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,J^k, (1) 



where Tj k is the distance from points r 3 - to r k , fj k is the unit vector from in 
the direction from points rj to r k , and 



A j:j = a j 



(2) 



where ay is the polarizability [15| of each dipole. Equation ([T]) represents the 
coupling between the dipoles and (T5]) allows the system of equations to be 
written in the compact form 



N 



^2 A i fcP i ~ E 



(3) 



since the dipole moment due to the total field at each dipole is P 3 = ocjEij. 
Because Aj k is a square matrix, the dipole moments can rapidly solved using 
the generalized minimum residual method 16]. 



3. Optimizing the Discrete Dipole Approximation interaction ma- 



The size of the microdevices we model may exceed 10-20 wavelengths in 
size, which may well require computational time in excess of several days 
and RAM beyond that available. To circumvent these limitations, we exploit 
the discrete rotational and/or mirror symmetry of a micro component. This 
is closely tied with the link between DDA and the T-matrix method. In 
the T-matrix method, the fields are represented as sums of vector spherical 
wavefunctions (VSWFs) p, and to use DDA to calculate a T-matrix, 
we can simply calculate the scattered field (and its VSWF representation) 
for each possible incident single-mode VSWF field in turn. The important 
point is that each VSWF is characterized by a simple azimuthal dependence 
of exp(im</>), where m is the azimuthal mode index. If we consider a group of 



Figure 2: Dipoles in a 4- fold discrete rotationally symmetric arrangement about beam 
axis. 

dipoles that are rotationally symmetric about the vertical axis — in the case 
of figure [U there is 4th-order rotational symmetry — the magnitude of the 
incident field will be the same, the field differing by the vector rotation Q 
and phase factor of exp(im0); only the dipole moment of one repeating unit 
of the total number of dipoles needs to be known. In spherical coordinates, 



trix 





(4) 
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However, our implementation of the interaction matrix[3], electric field and 
dipole moments were in cartesian coordinates system. This requires a trans- 
formation of the dipole moment at the first segment from the cartesian to 
spherical coordinate system followed by a transformation back from spherical 
to cartesian at the rotational counterpart dipole, 



\cart 



C g SP c r t exp(im < ; 
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cos(fl) 
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(6) 



is an orthogonal matrix that transforms the P^ art vector into the spherical 
coordinate system and 



sin(#) cos(0q) sin(6') sin( 
cos(#) cos(0q) cos(6 l ) sin( 
- sinO„) cos(d) n ) 



cos(fl) 
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is a transpose of S but at the azimuthal coordinate of the rotational coun- 
terpart dipole. The same transformations and phase corrections are applied 
to the E and H fields. Conventionally, the interaction matrix as defined in 





Figure 3: (a) 8-dipole example, (b) The 2 dipoles required to completely specify all dipole 
moments. 

[3] represents the coupling between each dipole with all other dipoles. Fig- 
ure HUa) shows the interaction matrix for the example set of dipoles shown 
in figure [3](a). In general the matrix will be made up of N x N cells for 
N dipoles; each cell is a 3 x 3 tensor. In the example, the matrix is made 
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up of 8 x 8 cells. A diagonal cell represents the self interaction (or the in- 
verse of the polarizability) and an off-diagonal cell represents the coupling 
between different dipoles. Taking advantage of the equal amplitudes and 
known phase factors between a dipole and its rotational counterparts, we 
can reduce the interaction matrix. Taking the example in figure [3]^b), we 
contruct the interaction matrix as if there were only 2 dipoles but we ag- 
gregate the contribution from the appropriate dipoles. For the off-diagonal 
cells, the coupling between a dipole with the other dipoles including their 
rotational counterparts are summed as follows 
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9=1 
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exp(i/cn 



' jk 



2(9) 
f jk 



x CgSexp(im0), j ^ k, 
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where Q is the order of discrete rotational symmetry, m is the azimuthal 
mode of the incident VSWF field, q is the rotational segment number, the 
rotational angle = 2-rrq/Q, rj^ is the distance from points rj to the ro- 

tationally symmetric points r k , and is the unit vector from points rj 

to The coordinate for a given rotational symmetric point is calculated 
using a rotation about the z-axis, 



,(«) 
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For the diagonal cells, the "self interaction" includes the coupling between a 
dipole and its rotational counterparts: 
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Figure IU(b) shows the interaction matrix representation for the example 
dipole system in figure El^c). The compressed interaction matrix is a fac- 
tor of Q 2 smaller than the conventional matrix. Having precalculated the 
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Figure 4: (a) Full interaction matrix, (b) Compressed interaction matrix. 



incident fields Ej iriC at each dipole of the rotational unit, we solve for the 
dipole moments Pj for the dipoles with a reduced set of linear equations: 



N 



k=l 



quad,j 



E 



inc,quad,j • 



(11) 



Notice that the discrete rotationally optimized A- matrix is m dependent; the 
A-matrices can be precalculated and loaded as cycle the required azimuthal 
modes (m) of the incident beam. However, there is some degree of redun- 
dancy; e.g., when m is even, A(m) = A(— m). The dipole moments of the 
rotational counterpart dipoles can be calculated by applying the rotational 
matrix [9] and phase factor exp(im</>). 

We can exploit mirror symmetry in a similar fashion, since a VSWF 
possess either even or odd parity with respect to the the xj/-plane. In the 
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Cartesian coordinate system, when n + m is odd, 
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and when n + m is even, 
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This allows a reduction in size of the A- matrix by a further factor of 4. 



4. The incident beam 

The optical tweezers that are used to drive micromachines are usually 
tightly focused laser beams. The beams my be linearly or circularly polarized 
(has spin angular momentum), and may or may not carry orbital angular 
momentum. For example, we use an LG 2 beam, which carries orbital angular 
momentum of 2h per photon, to trap and rotate the microrotors such as the 
one in figure [Jji. 

We calculate the incident field at each dipole using the vector spherical 
wave function expansion (VSWF) 



E - = EE a nmM®(kr) + b nm N®(kr), (14) 



n=l m= 



where k is the wave vector, r is the dipole position in spherical coordinates, 
n is the radial mode index, m is the azimuthal mode index and Mi 3 ! & Nl 3 i, 



are regular VSWFs [18|, [19(; a nm and b nm are incident coefficients for the 
illuminating beam calculated using incident beam functions from jij which 
uses the point matching method. If it sufficient to terminate the multipole 



expansion at n — N max where N max = ka + 3 v ka 20 



5. Formulating the T-matrix with point matching 

Here, we present the point matc hing method as an alternative to an ex- 
isting T-matrix method with DDA 13]. The point matching method that 



we implemented involved matching the fields due to the dipoles with the 
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fields calculated from vector spherical wave function expansion and at mul- 
tiple points. The number of points should be such that the linear system of 
equations are exactly or over determined i.e. the number of equations are 
the same or greater than the number of unknowns. 

Once we have the dipole moments, the E-field at any point r (relative 
to the origin) can be calculated by adding up the contributions from each 
dipole using the electric dipole field equation from section 9.2 of [2l| . 
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In the far field, the 1/r 2 and 1/r 3 terms can be ignored as their contributions 
dimish rapidly with distance: 
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5.1. Near field point matching 

The near field calculated via DDA f|T5|) can be matched to scattered field 
f l22|) calculated via the VSWFs at points around the scatterer (figure [5]). 
Since we have the dipole moments, Pj, we can derive a 'field matrix', F^-, 
derived from (fT5"|) such that 
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(DDA) 



Npz 



5> 



(17) 



where Npm is the number of point to be matched, j is the index for the 
dipole, i is the index for the near field point and each element Fij is a 3 x 3 
tensor. Using the vector identity (a x b) x c = b(a ■ c) — c(a • b), the cross 
product terms in ffTBI) become 



(f x p) x f 



p(r • r) — r{r ■ pj 
p-r(rp). 



(18) 



Substituting the result of ffTB"]) into f|T5|) we obtain 
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Figure 5: (a) Near field matching and (b) octant near field matching for a dipole model 
of a cube. 



which is can also be used to formulate the A-matrix ([T]). Thus the field 
matrix in (TPTD is 



_ exp(ifcrj fc ) 



k 2 (f jk f jk - 1 3 ) + -(3f jk r jk 



(20) 



where fj k is the unit vector between a dipole and a matched point. F^- can 
be precalculated and used when required. 

Rotational (jSJ) and mirror ( IT3T) symmetry optimizations may also be ap- 
plied to (ITT)) such that we only need to calculate the fields of the matched 
points in one octant: 

N PM /8 

-g(DDA,oct) _ p(oct)p(oct) ^21) 

i=l 

where only an the fields for an octant are calculated (although all the dipole 
moments are used). Thus reduces the number of equations by a factor of 64 
compare to (ITT)) . Now, since the VSWF expansion of the scattered field is 

00 n 

E sca = E P- M - W + ?n m N« (fcr). (22) 

n=l m=— ?i 

we can solve for the scattering coefficients 

p nm = M^l(kr)/B^l (23) 

10 



n - tvtW (kr\lVS JJ1JA > 

Hum — L ^nm\tvl ) / ^TM,nrn- 

As we cycle through each combination of n and m, we obtain the solutions 
for the scattering coefficients p nm and q nm which represent coupling between 
the n and m incident and scattered modes; p nm and q nm together make up 
one column of the T-matrix at a time. 

5.2. Far field matching 

In this method we match the VSWF expansion of the scattered field (122!) 
with that calculated via DDA in far field. We calculate the DDA field in a 
way similar to (j!7p but instead, calculating the far field 

Npm 

krEf DA) = £ F * P i> ( 24 ) 

i=l 

where field matrix, Fy, from (120]) is multiplied by kr and with the latter 
term omitted to give 

Fij = exp(ikr jk )k 3 (r jk f jk - 1 3 ). (25) 
Now the VSWFs in the expansion of the scattered field (122]) are 

M« (At) = N^l^M), (26) 



fcriV n 

iVn ( h^Ukr) - ) B nm (0, 0) (27) 



where iV n = 1/ y^n(n + 1), fo( 1 ' 2 ) are spherical Hankel functions and B nm , 
C nm & P nm are vector spherical harmonics: 

B nm (0,<f>) =rV 7 ;T(M) 

= vxc;(M) 

a // >> 7777 

^V^^^*^ 1 (28) 



ii 



C nm (OA) =Vx(ry™M) 



A 7777 >> f/ 

O^fiM*) ~^n(0,cf>), (29) 



P„ m (M)=*7n(M) (30) 

where 7™(#, 0) is the normalised scalar spherical harmonics. and N n m 

are the TE and TM multipole fields respectively. In the farfield limit 22|, |19 
the VSWFs 02SD & ([27D become 

Mffitr) = exp(iA;r)C nro (0, 0) (31) 

NSi(Ar) = ^(-i)" +1 exp(ifcr)B nm (^,0). (32) 
At 

The scattered can be further simplified since lim exp(ifcr) — > 1, 

1 ►00 

MaW^R^tM) (33) 
N«(fcr) = ^(-ir +1 B nm (^0), (34) 
Now the far field can also be calculated using 

/crE(0, 0) ^ ^ b nm B nm -\- C nm C nm -\- dnmP nm.) (35) 
nm 

where the coefficients b nm & c nm can be calculated (trapezoidal integration) 
using 

/ krE ■ B nm da 

Onm - r R , , , , {^) 

_ j krE ■ C nm da 

Cnm — r- p ^ , , { & < ) 

J '—'rim ' ^nm' 1 " 

where /crE is calculated via DDA; this is where the point matching occurs. 
The last term in fl35l) can be ignored because the radial component fl30|) is 
negligible in the far field. If we then match the terms in ( l22l . ( l35|) & (J33j) 
the scattering coefficients can be determined: 
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So, as with the near field matching, we cycle through the n and m modes 
and insert p nm vertically concatenated with q nm into the appropriate column 
of the T-matrix for the given mode. 



6. Mode redundancy and the order of discrete rotational symmetry 

We can also exploit the azimuthal dependence of exp(im0) of each indi- 
vidual VSWF. The discrete rotational symmetry effectively provides a pe- 
riodic boundary condition, determining the periodicity with respect to the 
azimuthal angle 0. According Floquet's theorem 0, in the case of the plane 



wave diffracting from a periodic grating, the wave vectors of the scattered 
modes are 

k S cat = knc + W, i = ... - 2, -1, 0, 1, 2, ... (40) 

where p is the order of discrete rotational symmetry of the scatterer. The 
spherical wave equivalent is such that the azimuthal modes of the scattered 
light that couple to a given incident mode is 

m scat = m inc + ip, i = ... - 2, -1, 0, 1, 2, .... (41) 

For a structure with four-fold discrete rotational symmetry, p = 4, thus 
m sca = mi nc , mi nc ± 4,mj nc ± 8, .... The number of unknowns in the linear 
system is reduced by a factor of 4. The calculation time saving for this 
scheme is about 2 to 3 orders of magnitude. 



7. Mode-reduced T-matrix 

This scheme may be deemed in breach of the true purpose of the T-matrix 
i.e. it should be good for any illumination at the wavelength for which it was 
formulated. However, if we know that the only beam that we are going to 
use is, say the LGq 2 and that the scatterer e.g. figure [T^l only spins along 
the beam axis, and stays on the beam axis, then the only relevant scattering 
modes when calculating the T-matrix would be m = 1,3. 

Bearing in mind that we 'move' the scatterer along the beam axis to 
find the equilibrium position by calculating the axial force jij; we perform 
repeated calculations with the mode-reduced T-matrix without sacrificing 
any accuracy. Formulating the mode-reduced T-matrix, in most cases, takes 
a fraction of the time taken for the full-mode T-matrix. 
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8. Results 
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Figure 6: a)Linear and b)log-log plots of the time required to assemble the interaction 
matrices for versus the number of dipoles (or equivalent for the 4-fold discrete rotational 
symmetric methodology 



8.1. Symmetry optimization 

The rotational symmetry optimized DDA scheme produces results that 
are no different (well below round-off errors) to that of the standard DDA 
implementation; there is no appreciable difference that would be visible when 
plotting comparison data for say the phase function or the extinction. In- 
stead, we present comparison test results for the calculation times for con- 
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Figure 7: a)Linear and b)log-log plots of the time required to solve for the dipole moments 
versus the number of dipoles 



structing the A-matrix and for solving the linear equations to obtain the 
dipole moments. 

Note that the full equivalent number of dipoles are plotted (figures [6] & [7]) 
for the case of the rotational symmetric algorithm and not just the number 
of dipoles in the quadrant. 

Although the size of the 4-fold rotationally symmetric interaction matrices 
are l/16th that of the conventional interaction matrix, the times required to 
assemble are comparable (figure [6]) and even slightly more for the former be- 
cause the algorithms (JT|) &; (|2]) still requires each component of the rotational 
counterpart in each quadrant to be calculated, including the transformation 
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Figure 8: The phase function plots on the left column are for the plane perpendicular to 
the (E-field) polarization of the incident plane wave whereas those on the right column 
are for the parallel plane. The plots for varying size parameters (2.5, 3, 3.5 & 4) are from 
top to bottom. The relative refractive index is 1.33. 
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Figure 9: The phase function for cubes with the refractive index of 1.5 and widths of 
a) 0.75A b) 1.25A. 
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matrices ([7J & ED for changing coordinate systems back and forth. Moreover, 
the rotationally symmetric A-matrices are azimuthal mode, m, dependent 
and thus requires 2N max + 1 number of compressed T-matrices; this number 
can reduced by exploiting A(m) = A(— m) for m even. 

The important advantage is that the memory footprint of the 4-fold dis- 
crete rotational symmetry optimized interaction matrix is l/16th that of its 
conventional counterpart. For the full interaction matrix, we encountered 
the problem of disk swapping on a PC with 4Gb of RAM just with > 5000 
dipoles and as can be seen on figure El we stopped short at < 5000 dipoles 
because the time taken increased sharply as available memory became scarce. 
By contrast, the symmetry optimized interaction matrix did not encounter 
such problems for the range of number of dipoles we tested with. 

The reduction (l/16th) of the number of linear equations when calcu- 
lating the dipole moments results very significant time savings (figure [7]). 
The comparison of the respective scaling with the number dipoles can be 
determined from the gradients of the log-log plot in figure [7b. 

8.2. Phase functions for the spheres and cubes 

To test the integrity of the optimization methods and the T-matrix via 
point matching, we compare against the results of the GLMT functions for 



the sphere [24j and cube [25|. The phase functions in figure [8] were generated 
from the scatterring coefficients obtained from the T-matrices; the results 
show good agreement against the Mie solution, allowing for stair casing errors 
when trying to approximate a sphere with a cubic lattice arrangement of 
dipoles with finite lattice spacing. 

The rotational symmetry algorithm was tested for cubes, since they pos- 
sess 4-fold rotational symmetry and analytical results can be calculated. The 
T-matrix via point matching was also applied and the results (figure E l) co m- 
pared favourably against those produced from the GLMT method in 
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8.3. Cross rotor torque calculations 

All the optimization methods and the point matching method for for- 
mulating the T-matrix were applied to calculating the forces and torques 
imparted on the microrotor (figure COt) by a tightly focuused and trapping 
LGq2 laser beam. We were able to determine the axial equilibrium position 
(figure EK) i-e. the point at which the axial force curve crosses the z-axis 
with negative gradient. The torque varies dramatically as we move along the 
beam axis but we are only interested in the torque at equilibrium. There are 
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a number of uncertainties in experimental measurement that we endeavour 
to narrow down but the results between the model and experiment agree 
within the error limits. 

9. Discussion 

The rotational and mirror symmetry algorithms do not introduce any er- 
rors (well below round-off errors) to the DDA calculation. The construction 
of the interaction matrix takes a large portion of the calculation time for 
the DDA method in general; in principle, the compressed interaction ma- 
trix offers no advantage here. Because the compressed interaction matrix 
is azimuthal mode dependent, it can be cumbersome when formulating the 
T-matrix; we load the precalculated matrices as required. However, there are 
two major advantages of the symmetry optimization scheme. Firstly, their 
memory footprint is significantly smaller; we are otherwise unable to model 
our microrrotors, given their size, on desktop PCs. Secondly, the time sav- 
ings when calculating the dipole moments more than make up for the extra 
time spent for the compressed interaction matrices. 

The mode redundancy algorithms give further time savings of about two 
orders of magnitude, as with the mode-reduced T-matrices where applicable. 

There are a number of factors that contribute to errors and uncertain- 
ties between the quantities of interest between experimental and modelling 
results. In our DDA model, we assume that the material is linear and non- 
absorbing. The lattice spacing is finite, which leads to stair casing errors and 
an inaccurate representation of the target. The characterics of the driving 
beam does invariably differ to some degree compare to that of the intended 
beam due to imperfections in the optical setup. Taking all these factors into 
consideration, the methods in described in this paper gives us sufficient accu- 
racy in modelling light scattering from the micromachines that we continue 
to design and develop. 
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